home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
LIBRARY
/
SHDK_2
/
SHCMPLX.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1992-05-13
|
13KB
|
485 lines
{$I SHDEFINE.INC}
{$I SHUNITSW.INC}
{$D-,L-}
unit ShCmplx;
{
ShCmplx
A Complex Arithmetic Unit
by
Bill Madison
W. G. Madison and Associates, Ltd.
13819 Shavano Downs
P.O. Box 780956
San Antonio, TX 78278-0956
(512)492-2777
CIS 73240,342
Copyright 1991 Madison & Associates
All Rights Reserved
This file may be used and distributed only in accord-
ance with the provisions described on the title page of
the accompanying documentation file
SKYHAWK.DOC
}
interface
uses
TPCRT,
TPDOS,
TPSTRING,
TpInline,
ShErrMsg;
type
{$IFNDEF Gen87}
extended = real;
{$ENDIF}
ComplexElement = extended;
ComplexBaseType = record
Re,
Im : ComplexElement;
end;
Complex = ^ComplexBaseType;
PFreeRec = ^TFreeRec;
TFreeRec = record
Next : PFreeRec;
Size : pointer;
end;
const
cxHaltErr : boolean = true; {Stop program execution on
non-I/O errors}
cxOK = 0;
cxDivByZero = 200;
cxBadMagnitude= 201;
procedure CmplxDispose(A : Complex);
{Replacement for the normal TP Dispose procedure. MUST be used to dispose
of variables of type Complex.}
procedure CmplxInit;
{Initializes the complex arithmetic package.}
procedure CmplxDeInit;
{De-initializes the complex arithmetic package, releasing all the ring
buffer heap space.}
function CmplxError : integer;
function Cmp2Str(A : Complex; Width, Places : byte) : string;
{Converts a complex value to a string of the form (Re + Im i)}
function CmpP2Str(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in radians.}
function CmpP2StrD(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in degrees.}
procedure CAbs(A : Complex; var Result : ComplexElement);
function CAbsF(A : Complex) : ComplexElement;
{Returns the absolute value of a complex number}
procedure CConj(A : Complex; Result : Complex);
function CConjF(A : Complex) : Complex;
{Returns the complex conjugate of a complex number}
procedure CAdd(A, B : Complex; Result : Complex);
function CAddF(A, B : Complex) : Complex;
{Returns the sum of two complex numbers A + B}
procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
function RxCF(A : ComplexElement; B : Complex) : Complex;
{Returns the complex product of a real and a complex}
procedure CSub(A, B : Complex; Result : Complex);
function CSubF(A, B : Complex) : Complex;
{Returns the difference between two complex numbers A - B}
procedure CMul(A, B : Complex; Result : Complex);
function CMulF(A, B : Complex) : Complex;
{Returns the product of two complex numbers A * B}
procedure CDiv(A, B : Complex; Result : Complex);
function CDivF(A, B : Complex) : Complex;
{Returns the quotient of two complex numbers A / B}
procedure C2P(A : Complex; Result : Complex);
function C2PF(A : Complex) : Complex;
{Transforms a complex in cartesian form into polar form}
{The magnitude will be stored in Result^.Re and the angle in Result^.Im}
procedure P2C(A : Complex; Result : Complex);
function P2CF(A : Complex) : Complex;
{Transforms a complex in polar form into cartesian form}
{The magnitude is stored in A^.Re and the angle in A^.Im}
procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
{Raises a complex (in polar form) to a real power}
{===========================================================}
implementation
VAR OUTPUT : TEXT;
const
MaxStack = 5;
StackTop = MaxStack - 1;
Pi = 3.14159265358979;
PiOver2 = Pi / 2.0;
TwoPi = Pi * 2.0;
F180OverPi = 180.0 / Pi;
SpConj : byte = 0;
SpSum : byte = 0;
SpDiff : byte = 0;
SpRProd : byte = 0;
SpProd : byte = 0;
SpQuot : byte = 0;
SpPolar : byte = 0;
SpCartsn : byte = 0;
SpPower : byte = 0;
CmpInited : boolean = false;
type
StackArray = array[0..StackTop] of Complex;
var
Conj,
Sum,
Diff,
RProd,
Prod,
Quot,
Polar,
Cartsn,
Power : StackArray;
CmpError : integer;
procedure CmplxDispose(A : Complex);
{Replacement for the normal TP Dispose procedure. MUST be used to dispose
of variables of type Complex.}
var
HP : pointer;
B : PFreeRec;
IsFree : boolean;
FRec : TFreeRec;
begin
HP := Normalized(HeapPtr);
A := Normalized(A);
B := Normalized(FreeList);
while (A <> Complex(B)) and (B <> HP) do begin
B := Normalized(B^.Next);
end;
if (A = Complex(B)) then begin
exit;
end;
Dispose(A);
end; {CmplxDispose}
procedure ChkInit;
begin
if not CmpInited then
RunErrorMsg(400, 'Unit ShCmplx not initialized.');
end;
procedure CmplxInit;
var
T1 : byte;
begin {CmplxInit}
for T1 := 0 to StackTop do begin
New(Cartsn[T1]);
New(Conj[T1]);
New(Diff[T1]);
New(Polar[T1]);
New(Power[T1]);
New(Prod[T1]);
New(Quot[T1]);
New(RProd[T1]);
New(Sum[T1]);
end;
CmpInited := true;
end; {CmplxInit}
procedure CmplxDeInit;
{De-initializes the complex arithmetic package, releasing all the ring
buffer heap space.}
var
T1 : byte;
begin {CmplxDeInit}
for T1 := StackTop downto 0 do begin
{Flush the heap space}
if Sum[T1] <> nil then
Dispose(Sum[T1]);
if RProd[T1] <> nil then
Dispose(RProd[T1]);
if Quot[T1] <> nil then
Dispose(Quot[T1]);
if Prod[T1] <> nil then
Dispose(Prod[T1]);
if Power[T1] <> nil then
Dispose(Power[T1]);
if Polar[T1] <> nil then
Dispose(Polar[T1]);
if Diff[T1] <> nil then
Dispose(Diff[T1]);
if Conj[T1] <> nil then
Dispose(Conj[T1]);
if Cartsn[T1] <> nil then
Dispose(Cartsn[T1]);
{Reset the heap pointers}
Conj[T1] := nil; Sum[T1] := nil; Diff[T1] := nil;
RProd[T1] := nil; Prod[T1] := nil; Quot[T1] := nil;
Polar[T1] := nil; Cartsn[T1] := nil; Power[T1] := nil;
end;
{Reset the ring buffer pointers}
SpConj := 0; SpSum := 0; SpDiff := 0;
SpRProd := 0; SpProd := 0; SpQuot := 0;
SpPolar := 0; SpCartsn := 0; SpPower := 0;
CmpInited := false;
end; {CmplxDeInit}
function CmplxError : integer;
begin
CmplxError := CmpError;
CmpError := cxOK;
end;
function Real2Str(A : ComplexElement; Width, Places : byte) : string;
var
T1 : string;
begin
Str(A:Width:Places, T1);
Real2Str := T1;
end;
function Cmp2Str(A : Complex; Width, Places : byte) : string;
begin
if A^.Im >= 0.0 then
Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' + ' +
Real2Str(A^.Im, Width, Places) + 'i)'
else
Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' - ' +
Real2Str(Abs(A^.Im), Width, Places) + 'i)';
end;
function CmpP2StrD(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in degrees.}
begin
CmpP2StrD := '('+Real2Str(A^.Re,Width,Places)+' at '+
Real2Str((A^.Im*F180OverPi),Width,Places)+#248')';
end;
function CmpP2Str(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in radians.}
begin
CmpP2Str := '('+Real2Str(A^.Re,Width,Places)+' at '+
Real2Str((A^.Im),Width,Places)+' rad)';
end;
procedure IncPtr(var StackName : StackArray; var StackPtr : byte);
begin
StackPtr := (StackPtr + 1) mod StackTop;
if StackName[StackPtr] <> nil then begin
Dispose(StackName[StackPtr]);
end;
New(StackName[StackPtr]);
end;
function CAbsF(A : Complex) : ComplexElement;
begin
CmpError := cxOK;
CAbsF := sqrt(A^.Re * A^.Re + A^.Im * A^.Im);
end;
procedure CAbs(A : Complex; var Result : ComplexElement);
begin
Result := CAbsF(A);
end;
function CConjF(A : Complex) : Complex;
begin
CmpError := cxOK;
ChkInit;
Conj[SpConj]^.Re := A^.Re;
Conj[SpConj]^.Im := -A^.Im;
CConjF := Conj[SpConj];
IncPtr(Conj, SpConj);
end;
procedure CConj(A : Complex; Result : Complex);
begin
Result^ := CConjF(A)^;
end;
function CAddF(A, B : Complex) : Complex;
begin
CmpError := cxOK;
ChkInit;
ChkInit;
Sum[SpSum]^.Re := A^.Re + B^.Re;
Sum[SpSum]^.Im := A^.Im + B^.Im;
CAddF := Sum[SpSum];
IncPtr(Sum, SpSum);
end;
procedure CAdd(A, B : Complex; Result : Complex);
begin
Result^ := CAddF(A, B)^;
end;
function RxCF(A : ComplexElement; B : Complex) : Complex;
begin
CmpError := cxOK;
ChkInit;
RProd[SpRProd]^.Re := A * B^.Re;
RPRod[SpRProd]^.Im := A * B^.Im;
RxCF := RProd[SpRProd];
IncPtr(RProd, SpRProd);
end;
procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
begin
Result^ := RxCF(A, B)^;
end;
function CSubF(A, B : Complex) : Complex;
begin
CmpError := cxOK;
ChkInit;
Diff[SpDiff]^ := CAddF(A, RxCF(-1.0,B))^;
CSubF := Diff[SpDiff];
IncPtr(Diff, SpDiff);
end;
procedure CSub(A, B : Complex; Result : Complex);
begin
Result^ := CSubF(A, B)^;
end;
function CMulF(A, B : Complex) : Complex;
begin
CmpError := cxOK;
ChkInit;
Prod[SpProd]^.Re := A^.Re * B^.Re - A^.Im * B^.Im;
Prod[SpProd]^.Im := A^.Im * B^.Re + A^.Re * B^.Im;
CMulF := Prod[SpProd];
IncPtr(Prod, SpProd);
end;
procedure CMul(A, B : Complex; Result : Complex);
begin
Result^ := CMulF(A, B)^;
end;
function CDivF(A, B : Complex) : Complex;
var
Denom : ComplexElement;
begin
CmpError := cxOK;
Denom := B^.Re * B^.Re + B^.Im * B^.Im;
if Denom = 0.0 then begin
CmpError := cxDivByZero;
ChkInit;
Quot[SpQuot]^.Re := 0.0;
Quot[SpQuot]^.Im := 0.0;
CDivF := Quot[SpQuot];
IncPtr(Quot, SpQuot);
exit;
end;
Quot[SpQuot]^ := RxCF(1.0 / Denom, CMulF(A, CConjF(B)))^;
CDivF := Quot[SpQuot];
IncPtr(Quot, SpQuot);
end;
procedure CDiv(A, B : Complex; Result : Complex);
begin
Result^ := CDivF(A, B)^;
end;
procedure C2P(A : Complex; Result : Complex);
{Transforms a complex in cartesian form into polar form}
{The magnitude will be stored in Result^.Re and the angle in Result^.Im}
begin
CmpError := cxOK;
Result^.Re := CAbsF(A);
if abs(A^.Re) > abs(A^.Im) then begin
{Use the tangent formulation}
Result^.Im := ArcTan(abs(A^.Im) / abs(A^.Re));
end
else begin
{Use the cotangent formulation}
Result^.Im := PiOver2 - ArcTan(abs(A^.Re) / abs(A^.Im));
end;
if A^.Im > 0 then begin {first or second quadrant}
if A^.Re < 0.0 then {Second quadrant}
Result^.Im := Pi - abs(Result^.Im);
end
else begin {third or fourth quadrant}
if A^.Re > 0.0 then {Fourth quadrant}
Result^.Im := TwoPi - abs(Result^.Im)
else begin {Third Quadrant}
Result^.Im := Pi + abs(Result^.Im);
end;
end;
if Result^.Im = TwoPi then Result^.Im := 0.0;
end;
function C2PF(A : Complex) : Complex;
begin
ChkInit;
C2P(A, Polar[SpPolar]);
C2PF := Polar[SpPolar];
IncPtr(Polar, SpPolar);
end;
procedure P2C(A : Complex; Result : Complex);
{Transforms a complex in polar form into cartesian form}
{The magnitude is stored in A^.Re and the angle in A^.Im}
begin
CmpError := cxOK;
Result^.Re := A^.Re * cos(A^.Im);
Result^.Im := A^.Re * sin(A^.Im);
end;
function P2CF(A : Complex) : Complex;
begin
ChkInit;
P2C(A, Cartsn[SpCartsn]);
P2CF := Cartsn[SpCartsn];
IncPtr(Cartsn, SpCartsn);
end;
function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
{Raises a complex (in polar form) to a real power, returning the result
also in polar form}
begin
CmpError := cxOK;
if A^.Re <= 0.0 then begin
CmpError := cxBadMagnitude;
ChkInit;
Power[SpPower]^.Re := 0.0;
Power[SpPower]^.Im := 0.0;
exit;
end;
Power[SpPower]^.Re := exp(R * ln(A^.Re));
Power[SpPower]^.Im := R * A^.Im;
CpPwrRF := Power[SpPower];
IncPtr(Power, SpPower);
end;
procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
begin
Result^ := CpPwrRF(A, R)^;
end;
end.